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Abstract 

The swimming of an assembly of rigid spheres immersed in a viscous fluid of infinite extent 
is studied in low Reynolds number hydrodynamics. The instantaneous swimming velocity and 
rate of dissipation are expressed in terms of the time-dependent displacements of sphere centers 
about their collective motion. For small amplitude swimming with periodically oscillating dis¬ 
placements, optimization of the mean swimming speed at given mean power leads to an eigenvalue 
problem involving a velocity matrix and a power matrix. The corresponding optimal stroke per¬ 
mits generalization to large amplitude motion in a model of spheres with harmonic interactions and 
corresponding actuating forces. The method allows straightforward calculation of the swimming 
performance of structures modeled as assemblies of interacting rigid spheres. A model of three 
collinear spheres with motion along the common axis is studied as an example. 
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1. INTRODUCTION 


In earlier work jl[ we presented a method to analyze the performance of a microswim¬ 
mer modeled as an assembly of N rigid spheres immersed in a viscous incompressible fluid 
of inhnite extent, with a no-slip boundary condition on the surface of each sphere. The 
motion of the whole system is determined by the Stokes equations of low Reynolds number 
hydrodynamics. The swimming motion of such a system was discussed earlier by Alouges 
et ah The particular case of collinear spheres was studied by Vladimirov ^ using a 

two-timing method. 

For small displacements of the spheres from fixed positions in the collective rest frame 
the time-averaged swimming velocity and rate of dissipation can be evaluated in terms of a 
{3N — 3) X {3N — 3) velocity matrix and a {3N — 3) x {3N — 3) power matrix, which can be 
constructed from the mobility matrix for each relative rest configuration [I|. Optimization 
of the velocity at fixed power leads to a generalized eigenvalue problem involving the two 
matrices. Optimal efficiency corresponds to the maximum eigenvalue. 

In a model with harmonically interacting spheres the optimal stroke of small amplitude 
motion can be used to calculate a set of corresponding actuating forces. Large amplitude 
motion can be studied by solving the equations of Stokesian dynamics for the same actuating 
forces multiplied by a factor. The mean swimming velocity and the mean power of the large 
amplitude motion can then be determined numerically from the limit cycle of the solution. 

In the following we present an alternative method based on a purely kinematic point of 
view. Expressions are derived for the instantaneous swimming velocity and power in terms 
of the sphere displacements from the center and their instantaneous time derivative. This 
allows calculation of the mean swimming velocity and mean power for given periodic stroke 
of any amplitude. The present method also provides an alternative derivation of the velocity 
matrix and power matrix of small amplitude motion. 

For large amplitude swimming the present method is more straightforward than the earlier 
one [I| , since it does not require numerical solution of the equations of Stokesian dynamics. 
A large amplitude stroke may be determined by amplifying the optimal stroke found from 
the eigenvalue problem of the small amplitude theory for a given equilibrium structure. The 
instantaneous swimming velocity and power are then determined from explicit expressions 
in terms of the given displacements. Subsequently the mean swimming velocity and mean 
power can be found by integration over a period. 

Both methods are tested on a model of three collinear spheres with motion along the 
common axis, as formulated by Najah and Golestanian and studied in detail by Golesta- 
nian and Ajdari [^. The two methods of calculation lead to similar numerical results for a 
wide range of amplitude. 


2. DISPLACEMENT AND SWIMMING VELOCITY 

We consider a set of N rigid spheres of radii oi,..., oat immersed in a viscous incompressible 
fluid of shear viscosity r]. The fluid is of inhnite extent in all directions. At low Reynolds 
number and on a slow time scale the how velocity v and the pressure p satisfy the Stokes 
equations 

rj'V'^v — Vp = 0, V ■ r; = 0. (2.1) 
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The flow velocity v is assumed to satisfy the no-slip boundary condition on the surface of 
the spheres. The fluid is set in motion by time-dependent motions of the spheres. At each 
time t the velocity held v{r,t) tends to zero at infinity, and the pressure p{r,t) tends to 
the constant ambient pressure po- We shall study periodic relative motions which lead to 
swimming motion of the collection of spheres. 

We assume that the motion is caused by time-dependent periodic forces Fi(t),FN(t) 
which satisfy the condition that their sum vanishes at any time. The forces are transmitted 
by the spheres to the fluid. The spheres can rotate freely, so that they exert no torques on 
the fluid. Hence the rotational velocities f2i(t),..., f2Ar(t) can be ignored. The translational 
velocities Ui, ...,Un are linearly related to the forces, 

N 

J = (2.2) 

k=l 


with translational mobility tensors ^*4- The tensors have many-body character and depend 


in principle on the positions of all particles [8[-[10|- By translational invariance only relative 


distance vectors {Ri — Rj} occur in the functional dependence. We abbreviate eq. (2.2) as 

U = /X ■ F, (2.3) 


with a symmetric 3N x 3N mobility matrix /j,. Conversely 


F = C • U, (2.4) 

with friction matrix The friction matrix is the inverse of the mobility matrix, ^ 
and is also symmetric. 

The positions of the centers change as a function of time. The equations of motion of 
Stokesian dynamics read 

dR 

= C/,(i^l, ..., Rn, t), J = 1,..., N. (2.5) 

The explicit time-dependence on the right originates in the time-dependence of the forces 
F(t). In the swimming motion the forces are periodic in time with period T, so that F(t-|-T) = 
F(t). As mentioned, we impose the condition that at no time is there a net force acting on 
the set of spheres, so that 

N 

^F,(i) = 0. (2.6) 

i=i 

We look for a solution of eq. (2.5) corresponding to swimming motion, of the form 

R,{t) = Sjo+[ U{t')dt' + 6j{t), j = l,...,iV, (2.7) 

Jo 

where the first two terms describe the collective motion of the configuration So = 
{Sio,Sp^o) with swimming velocity U(t) caused by the displacements {(5j(t)}. We re¬ 
quire that the latter are periodic with period T, and exclude uniform displacements, so that 
the 377-dimensional vector d(t) = {5i(t),..., Sp[(t)} satisfies 

d{t) ■ u„ = 0, (a = X, p, z), (2.8) 
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where the symbol Ux denotes a 3iV-dimensional vector with 1 on the x positions, 0 on the 
y, z positions, and cyclic. Periodicity implies 

U{t + T) = U{t), d(t + T) = d(t). (2.9) 


The mean swimming velocity is dehned as 


- 1 r 

Usw = 7^ J u{t) dt. 

We require that d(t) is purely oscillating, so that 



0 . 


( 2 . 10 ) 


( 2 . 11 ) 


We show in the following that the instantaneous swimming velocity U{t) can be calcu¬ 
lated from the displacement vector d(t) and its time derivative d(t). Later we compare the 
present kinematic description to a dynamical model, in which the forces are decomposed 
into actuating forces and elastic restoring forces. 


3. SWIMMING VELOCITY AND DISSIPATION 


By substitution of eq. (2.7) into eqs. (2.4) and (2.5) one finds 

F = C-(f//3U;3 + d), 


(3.1) 


where summation over repeated greek indices is implied. The condition (2.6) can be ex¬ 
pressed as Uq, ■ F = 0, so that 

ZayUp = —Ua ■ C ■ d (3.2) 

with friction tensor 

Zap LIq. ■ ^ ■ Uy. (3.3) 

Hence we obtain the swimming velocity 


Uci ^ ■ d. 


(3.4) 


where May is the inverse of the friction tensor. The 3N x 3N friction matrix ^ depends only 
on the instantaneous relative positions. Therefore the friction tensor Z and the mobility 
tensor M depend on the displacement vector d, but not on the central coordinates Ra = 

■ R/N. 

By series expansion of the mobility tensor M and the friction matrix ^ in powers of the 
displacement vector d we obtain a corresponding expansion of the swimming velocity 

Lr = t/(i) + t/(2) + c/(3) + ..., (3.5) 

with hrst order term 

Clt’ = ■ C“ ■ d, (3.6) 

with mobility tensor and friction matrix ((('^ calculated for the conhguration Sq. By 
periodicity of d(f) the time average of the first order swimming velocity vanishes, = 0. 
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We introduce the friction vectors 


fa = u„ ■ C = C ■ Ua, (3.7) 

where we have used the symmetry of the friction matrix The vectors are related to the 
friction tensor by 

Ua ■ f/3 ■ fee Zee/S- (3.8) 

From the Taylor series expansion of eq. (3.4) we hnd that the second order instantaneous 
swimming velocity can be expressed as 



t/f = -d.v“yd, 

(3.9) 

with matrix V" given by 

V“ = 

(3.10) 

where V is the gradient vector in 3iV-dimensional configuration space. The notation in 
eq. (3.9) indicates that the matrix-function is to be evaluated at R = Sq. 

The expression on the right of eq. (3.10) may be written as a sum of two terms, 


V“ = (VM„^)f^ + M,^D^ 

(3.11) 

with derivative friction matrix 

= Vf^. 

(3.12) 

We introduce the gradient vectors 



g? = ■ u, = 

(3.13) 

and use the identity 


(3.14) 

to show that 

= —MaygjMs/s- 

(3.15) 

Then eq. (3.11) may be expressed alternatively as 




(3.16) 

with reduced derivative friction matrix 



6^ = - g^^M^sfs- 

(3.17) 

This matrix has the property 

6^ ■ uq, = 0. 

(3.18) 

From the fact that ^ depends only on relative coordinates it follows that Uq, ■ = 

hence 

Ua ■ = 0, Uce-g^ = 0. 

0, and 

(3.19) 

As a consequence 

o 

> 

o' 

> 

(3.20) 
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The time-dependent rate of dissipation can be expressed in the same matrix formalism. 
The rate of dissipation is given by 

P=F-U = F-cI, (3.21) 

since F ■ Uq, = 0 on account of the condition eq. (2.6). Substituting eq. (3.1) we hnd 

P = d ■ C ■ d + f/«cl ■ f„. (3.22) 

It follows from eq. (3.4) that the rate of dissipation is at least of second order in d and d. 
To second order, by use of eq. (3.6), 


= d ■ P . d (3.23) 

with matrix 

P = c° - M°^f°f°. (3.24) 

The matrix is symmetric and has the properties 

u„ ■ P = 0, P ■ u„ = 0. (3.25) 

The properties eq. (3.20) and (3.25) allow us to reduce the dimension of the matrix descrip¬ 
tion by three by the introduction of center and relative coordinates. 


4. VELOCITY MATRIX VECTOR AND POWER MATRIX 

The center of the assembly is given by 

1 ^ 1 

-R= ^ ■ R (4.1) 

i=i 

with Cartesian unit vectors e^. We dehne relative coordinates {r'j} as 

ri = i?2 — Rl, ^2 = i?3 — i?2, •••, 

'f‘N-1 = Rn — Rn- 1, j = — (4.2) 

and the corresponding (3iV — 3)-vector r = (ri,..., The 3iV-vector {R, r) is related to 

the vector R by a transformation matrix T according to 

{R,r) = J-R (4.3) 

with explicit form given by eqs. (4.1) and (4.2). 

The matrices V" and P are transformed to 

V“ = T-V“-T-\ Pt = T .p .T-\ (4.4) 

The first three rows of T consist of u^/N and the first three columns of T“^ consist of Uq. 
It follows from the properties eq. (3.20) and (3.25) that the hrst three rows and columns of 
the transformed matrices and Pr vanish identically. Hence in this representation we can 
drop the center coordinates and truncate the matrices by erasing the hrst three rows and 
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columns. We denote the truncated (3iV — 3) x (3iV — 3)-matrices as and Pr and define 
displacements ^ in relative space by 


(0,O=T-d. (4.5) 

With this notation the second order swimming velocity and rate of dissipation are given by 

t/f = I. Ct ■ ■ €, i . Ct ■ Pr ■ (4.6) 

with the matrix 

Ct= [T-i-T-V- (4.7) 

This (3iV —3) X (3iV —3) dimensional matrix consists of numerical coefficients and is obtained 

from the corresponding 3N x 3N matrix by truncation, as indicated by the final hat symbol. 
We consider in particular harmonically varying displacements of the form 

d(t) = ds sincat + dcCoscaf, (4.8) 

with a corresponding expression for ^{t). The time-averaged second order swimming velocity 


and rate of dissipation are then given by 

^ ( 4 . 9 ) 

We introduce the complex dimensionless vector 

+ ( 4 . 10 ) 

where 6 is a typical length scale. With the definitions 

= A=Fc^.Pt, (4.11) 

and the scalar product 

N-l 

(€V) = E «? ■ F (4.12) 

t=l 

the mean swimming velocity and mean rate of dissipation can then be expressed as 

W = ^cu 6 (riB“ir), ^ (4.13) 


We have normalized such that the matrix elements of B" and A are dimensionless. We call 
B“ the velocity matrix and A the power matrix. 

We ask for the stroke with maximum swimming velocity in a class of strokes with equal 
rate of dissipation for fixed values of the geometric parameters, fixed frequency u, and fixed 
viscosity 77 . This leads to the generalized eigenvalue problem 

= A“A|fi (4.14) 
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The eigenvalues {A"} are real. The maximum efficiency for motion in direction a is given 
by the maximum eigenvalue as 

= A”„. (4,15) 

The set -^rmaan ^Tmax} depends on the choice of Cartesian coordinate system. Fur¬ 

ther optimization may be possible by a rotation of axes. In particular cases a natural choice 
of axes will suggest itself. 

In the formulation of the mobility matrix in Eg. (2.2) the nature of the forces {E'j} 
need not be specihed. In an earlier calculation jll| we considered microswimmers with 
internal harmonic interactions, driven by actuating forces. In matrix form the forces may 
be expressed as 

F = E + H-(R-So), (4.16) 

with a real symmetric matrix H with the property H ■ Uq, = 0. The actuating forces 
are assumed to satisfy 

N 

(4,17) 

i=i 

They can be generated internally or externally. 


5. THREE-SPHERE SWIMMER 


The simplest application of the theory is to a three-sphere swimmer with three spheres 
aligned on the x axis, as studied by Golestanian and Ajdari [^. The spheres move along the 
X axis, and the y and z coordinates can be ignored. There are only two relative coordinates 
Ti = X 2 — Xi and r 2 = x^ — X 2 , and the relevant parts of the matrices and A are two- 
dimensional. The elements of the 3x3 mobility matrix are approximated by use of the 
Oseen interaction as 




6717] 


Sjk + 


2\x-j - Xk 


(1 — Sjk) 


(5.1) 


In the bilinear theory we consider a point ro in r-space with coordinates {di,d 2 ), corre¬ 
sponding to the conhguration Sq of the rest system. As an example we consider the case of 
equal-sized spheres with oi = 02 = 03 = a and equal distances between centers di = d 2 = d. 
For this case the explicit expressions for the matrices and A are identical to those derived 
earlier by a different method [l[ . Explicit expressions for the eigenvectors and eigenvalues 
A± of the two-dimensional eigenvalue problem B^ • ^ = AA^, as functions of the ratio d/a, 
were derived in ref. 1 . 

In the bilinear theory, corresponding to small £, the orbit (ri(f),r 2 (t)) = {x 2 {t) — 
Xi{t),X 3 (t) — X 2 (t)) in relative space is given by r(t) = ^ 0 - 1 - ^o(^) with Vq = (d, d) and 


^o(^) = Re exp{—iojt), 


(5,2) 


with amplitude factor e and eigenvector = ( 1 , ,^+) corresponding to the largest eigenvalue. 
In Eg. 1 of ref. 1 we have shown the elliptical orbit in relative space for d = 5a and e = 0.1. 
The corresponding displacement vector in configuration space is given by 


0 

Ut) 





Mt) = T-' ■ 


(5.3) 






In fig. 1 we show the reduced mean swimming velocity Usw /as a function of e for 
d = 5a, as calculated from eq. (3.4). In £g. 2 we show the reduced mean rate of dissipation 
V /as calculated from eq. (3.22). In fig. 3 we show the efficiency Et = r]ua‘^Usw/T^ 
as a function of e. The efficiency increases monotonically with the amplitude factor. 

It is of interest to compare the above results with values obtained by the numerical 
solution of the Stokesian equations of motion eq. (2.5) with hydrodynamic interactions 
given by eq. (5.1) and prescribed oscillating actuating forces. We use harmonic interactions 
given by the 3 X 3-matrix 

/-I 1 0 \ 

H = A; 1 -2 1 (5.4) 

VO 1-1/ 


with elastic constant k. This corresponds to nearest neighbor interactions of equal strength k 
between the three spheres. The stiffness of the swimmer is characterized by the dimensionless 
number a defined by 


(5.5) 


Tirjaoj 


In general, the hrst order forces Fq^V^) corresponding to the displacement vector do(t) and 
the corresponding first order swimming velocity calculated from eq. (3.6), follow 

from eq. (3.1) as 




(5.6) 


In the present case only the x components are relevant. The corresponding actuating forces 
Eo(t) are found from eq. (4.16) as 


£„(«) = F|,'>(«)-H-d„(«). 


(5.7) 


These have the property Uq, ■ Eo(t) = 0, so that the sum of actuating forces vanishes. We 
choose initial conditions for the x coordinates 


a;i(0) = 0, X2{0) = d + ea, 0 ) 3 ( 0 ) = 2(i-(-ea-|-ea Re (5.8) 

In fig. 4 we show the numerical solution of the equations of Stokesian dynamics eq. (2.5) 
with forces given by 

F(t) = Eo(t) + H-(R(f)-So) (5.9) 

for d = 5a, stiffness cr = 1, and amplitude factor e = 2 for the hrst ten periods. We compare 
the orbit with the ellipse given by eq. (5.2). The mean swimming velocity and mean power, 
calculated as time-averages over the last period for values of the amplitude factor in the 
range 0 < e < 2, are shown in kgs. 1 and 2. The corresponding efficiency is shown in 
kg. 3. The dashed curves in kgures 1 — 3 replace kgs. 3, 4, and 5 of ref. 1, which were 
calculated from inappropriate actuating forces. The efficiency is approximately twice as 
large as calculated in ref. 1 . 

It is true that in kg. 3 the efficiency for given e calculated by the kinematic method 
is always larger than that calculated by the dynamic method from the limit cycle with 
actuating forces. However, we must compare the mean swimming velocity for two dikerent 
strokes of the same mean power. In kg. 5 we plot the power as a function of e in the range 
1.9 < e < 2 as calculated by the two different methods. The value D = 52 T]uj‘^a^ of the 
mean power occurs at Sk = 1.949 in the kinematic method, and at Ed = 1.970 in the dynamic 
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method. For these values the mean swimming velocity is found to be Usw = 0.0546 ua for 
the elliptical orbit of the kinematic method, and Usw = 0.0538 uja for the limit cycle of the 
dynamical method. Thus in the present case the elliptical orbit is the most efficient of the 
two. This does not exclude that for the same power an orbit with yet higher speed can be 
found. 

At £ = 1.38 and for d = 5a we have Usw ~ 0.026 ua from eq. (3.4) and V ^ 25.8 rju'^a^ 
from eq. (3.22) for the orbit given by eq. (5.3). This can be compared with the numerical 
calculation of Alouges et ah [2|,j^ on the basis of a Stokes solver. The authors used radius 
a = 0.05 mm, and period T = 1 s. For viscosity of water r] = 0.01 poise our calculation 
yields A = UswT ~ 0.0081 mm and VT ^ 0.127 x 10“^^ J. The latter value is somewhat 
less than the one given in table 1 of ref. 3, and the displacement agrees well with the value 
0.01 mm of Alouges et ah. 

Finally we consider the efficiency calculated from eqs. (3.4) and (3.22) for displace¬ 
ment in relative space of the form eq. (5.2), but with the eigenvector replaced by 
^ = (1, Aexp(id)) with absolute value A and phase S. The values of A and S can be related 
to the Stokes parameters of the elliptical orbit [l^. In fig. 6 we show the efficiency for 
amplitude factor e = 2 and ratio d/a = 5 as a function of A and S. The maximum is not 
very pronounced. 


6. DISCUSSION 

The swimming performance of an assembly of spheres as a function of the amplitude of 
a chosen stroke can be studied in a purely kinematic formulation. From eq. (3.4) we find 
the instantaneous swimming velocity, and from eq. (3.22) we find the instantaneous rate of 
dissipation or power. The mean swimming velocity and the mean power follow by averaging 
over a period. The ratio of these two quantities yields the efficiency of the stroke. 

Alternatively one may use a dynamic approach [H,[lH in which the swimmer is modeled as 
a set of spheres bound harmonically to equilibrium positions and with harmonic interactions. 
The spheres are subject to actuating forces which sum to zero. The corresponding swimming 
motion may be found as the limit cycle of the solution of the equations of Stokesian dynamics. 
The mean swimming velocity and the mean power may be found numerically from the limit 
cycle. 

We have shown in sect. 5 that for a collinear three-sphere swimmer the two methods 
lead to similar results over a wide range of amplitude, provided that for small amplitude the 
actuating forces correspond to the chosen kinematic stroke. We have chosen the latter to be 
the optimal one at small amplitude, as determined from the velocity matrix and the power 
matrix of the bilinear theory. 

The kinematic method is the more straightforward one, since it does not require numerical 
solution of the equations of Stokesian dynamics. The dynamic approach has the advantage 
that it provides a physical model of the swimmer. It will be of interest to explore the 
difference in efficiency for given stroke or given actuating forces as a function of amplitude 
factor for more sophisticated model swimmers, with actuating forces chosen to agree with 
the optimal stroke at small amplitude. 
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Figure captions 


Fig. 1 

Plot of the reduced mean swimming velocity Uaw/[s^ojo) for d = 5a as a function of 
the amplitude £ as calculated by the kinematic method (solid curve), and by the dynamic 
method with stiffness parameter <7=1 (dashed curve). 

Fig. 2 

Plot of the reduced mean swimming power V /for d = 5a as a function of 
the amplitude e as calculated by the kinematic method (solid curve), and by the dynamic 
method with stiffness parameter <7 = 1 (dashed curve). 


Fig. 3 

Plot of the efficiency = r]ua^Usw/T^ for d = 5a as a function of the amplitude e as 
calculated by the kinematic method (solid curve), and by the dynamic method with stiffness 
parameter <7 = 1 (dashed curve). 


Fig. 4 

Plot of the orbit in the rir 2 plane calculated from the equations of Stokesian dynamics 
for d = 5a, e = 2, <7 = 1 for ten periods. The initial values correspond to Eq. (5.8) and 
the forces follow from eq. (5.9). We also plot the elliptical orbit for d = 5a, e = 2 (dashed 
curve). 


Fig. 5 

Plot of the mean swimming power V/{r]u‘^a^) for d = 5a as a function of the amplitude 
e in the range 1.9 < e < 2 as calculated by the kinematic method (solid curve), and by the 
dynamic method with stiffness parameter <7 = 1 (dashed curve). 


Fig. 6 

Plot of the efficiency Et = rjua^Usw/T^ calculated by the kinematic method for the 
elliptical orbit in the rir 2 plane given by eq. (5.2) for d = 5a with e = 2 and replaced 
by ^ = (1, Aexp(id)) as a function of amplitude A and phase 6. 
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